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ABSTRACT 

We fitted Spitzer/IRS ~ 2 — 35 flm spectra of 26 luminous QSOs in attempt to define the main emission 
components. Our model has three major components: a clumpy torus, dusty narrow line region (NLR) clouds 
and a blackbody-like dust. The models utilize the clumpy torus of Nenkova et al. (2008) and are the first to 
allow its consistent check in type-I AGNs. Single torus models and combined torus-NLR models fail to fit 
the spectra of most sources but three component models adequately fit the spectra of all sources. We present 
torus inclination, cloud distribution, covering factor and torus mass for all sources and compare them with 
bolometric luminosity, black hole mass and accretion rate. The torus covering factor and mass are found to be 
correlated with the bolometric luminosity of the sources. We find that a substantial amount of the ~ 2 — Ifim 
radiation originates from a hot dust component, which likely situated in the innermost part of the torus. The 
luminosity radiated by this component and its covering factor are comparable to those of the torus. We quantify 
the emission by the NLR clouds and estimate their distance from the center. The distances are ~ 700 times 
larger than the dust sublimation radius and the NLR covering factor is about 0.07. The total covering factor by 
all components is in good agreement with the known AGN type-I:type-II ratio. 
Subject headings: infrared: galaxies - galaxies: active - galaxies: nuclei - quasars: general 



1. INTRODUCTION 

An obscuring dusty structure surrounding an accreting mas- 
sive black hole (BH) is believed to be a common feature of 
most active galactic nuclei (AGNs). Since the obscuration of 
the central region is anisotropic, sources with small inclina- 
tion angles to the line of sight (face-on) would be classified 
as type-I AGNs and those with large inclination angles (edge- 
on) would be classified as type-II AGNs. In this picture, the 
bulk of the radiation from the central engine is absorbed by 
the obscuring structure and re-emitted mainly in mid-infrared 
(MIR) wavelengths. The central obscuration is not necessar- 
ily a single component structure and the exact nature of its 
different components remains an open question and provides 
the motivation to the present work. 

A main component of the obscuring structure is believed 
to be a dusty torus. The MIR spectral energy distribution 
(SED) of such torus depends on its dimensions and geometry, 
the density distribution and the dust grain properties. Initial 
attempts to model such tori assumed smooth density distri- 
butions (e.g. Pier & Krolik 1992; 1993; Granato & Danese 
1994; Efstathiou & Rowan-Robinson 1995; van Bemmel & 
Dullemond 2003; Schartmann et al. 2005). This can explain 
part of the SED but falls short of providing realistic MIR spec- 
tra. For example, Silicate emission has been detected in Type- 
2 AGNs (e.g. Sturm et al. 2006; Teplitz et al. 2006) despite 
of the the fact that absorption features are expected in edge-on 
tori. Several studies suggested that the dusty medium should 
be clumpy (e.g. Krolik & Begelman 1988; Rowan-Robinson 
1995; Nenkova et al. 2002; Tristram et al. 2007; Thomp- 
son et al. 2009). The recent works of Nenkova et al. (2002) 
and Nenkova et al. (2008a; 2008b; hereafter N08) offer the 
required formalism to calculate the SED of such tori and a 
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framework to constrain some of their properties. 

Other AGN components can contribute to the observed 
MIR spectrum of AGNs. Some of this emission may originate 
farther from the central radiation source, at distances exceed- 
ing the dimensions of the torus. Broad-band 10 ftm imaging 
of several nearby AGNs suggests extended MIR continuum 
source (Cameron et al. 1993; Bock et al. 2000; Tomono et 
al. 2001; Radomski et al. 2003; Packham et al. 2005). Dusty 
clouds in the narrow line region (NLR) may be the source of 
such radiation (Schweitzer et al. 2008; hereafter S08). Thus, 
a significant contribution at ^10-30 flm due to components 
not related to the torus, must be considered. 

Another component not necessarily related to the torus is 
hot dust emission in the immediate vicinity of the central en- 
gine. Minezaki et al. (2004) reported delayed and corre- 
lated variations between the V and the K band emission in 
the nucleus of the Seyfert 1 galaxy NGC4151. The mea- 
sured lag between the V and K bands, ~ 48 days, lead to 
the conclusion that the near infrared emission is dominated 
by thermal radiation from hot dust ^0.04 pc from the cen- 
ter. This result was recently supported by the work of Riffel 
et al. (2009) who fitted the near infrared spectrum of NGC 
4151 with ~ 1300/f blackbody (BB) spectrum representing 
emission from hot dust in the inner region of the torus. Other 
studies used similar models to fit the SED of more luminous 
AGNs (e.g. Edelson & Malkan 1986; Barvainis 1987; Kishi- 
moto et al. 2007). The study of Suganuma et al. (2006) indi- 
cates that luminous and variable K-band emission is common 
in several nearby Seyfert- 1 galaxies. More generally, obser- 
vational support for powerful 1-3 flm emission is known for 
decades (e.g., Hyland & Allen 1982; McAlary et al. 1983; 
Neugebauer et al. 1987). It is not clear that such emission 
is consistent with the torus dust emission, which is expected 
to peak at longer wavelengths. Furthermore, part of this radi- 
ation requires dust temperature that exceeds the sublimation 
temperature of silicate type grains suspected to be responsible 
for the bulk of torus MIR radiation. 
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TABLE 1 
The QUEST QSO sample 



Object 


z 


Dl 


log Lgioo 


NIR Data 








(PYP S ^ t 


References 


PG 2349-014 


0.1740 


840 


44.81 


1,2,3,4,5,6 


PG 2251+113 


0.3255 


1706 


45.63 


3,4,7,8 


PG 2214+139 


0.0658 


295 


44.40 


2,4,6,7,9,10,11,12 


PG 1700+518 


0.2920 


1504 


45.68 


3,4 


PG 1626+554 


0.1330 


623 


44.44 


4 


PG 1617+175 


0.1124 


520 


44.29 


3,4 


PG 1613+658 


0.1290 


603 


44.70 


1,2,3,4,7,11,13 


PG 1448+273 


0.0650 


291 


44.283 


2,3 


PG 1440+356 


0.0791 


357 


44.22 


3,4,10,13 


PG 1435-067 


0.1260 


589 


44.39 


7 


PG 1426+015 


0.0865 


393 


44.44 


2,3,4,7 


PG 1411+442 


0.0896 


408 


44.31 


3,4,13 


PG 1309+355 


0.1840 


891 


45.081 


4 


PG 1302-102 


0.2784 


1424 


45.17 


4,14 


PG 1244+026 


0.0482 


213 


43.593 


2,7 


PG 1229+204 


0.0630 


28 1 


43.895 


1,2,3,4,13,15,16 


PG 1126-041 


0.0600 


267 


43.82 


1,2,4,6,13,16 


PG 1116+215 


0.1765 


851 


45.13 


4,6,7,17 


PG 1004+130 


0.2400 


1201 


45.42 


3,4,8,14,18 


PG 1001+054 


0.1605 


766 


44.55 


4,8,13,14 


PG 0953+414 


0.2341 


1168 


45.11 


3,4 


PG 0838+770 


0.1310 


613 


44.16 


2,3,4,13 


PG 0026+129 


0.1420 


670 


44.66 


3,4,7,8,14,18 


PG 0157+001 (Mrk 1014) 


0.1630 


779 


44.68 


2,3,4,5,16,19,20,21 


PG 0050+124 (IZwl) 


0.0611 


273 


44.30 


6,7,9,10,12,16,22,23,24 


B2 2201+31 A 


0.295 


1522 


45.91 


4,8 



REFERENCES. — (1) Rudy et al. (1982); (2) 2MASS magnitudes (Jarrett et al. 2000); (3) Neugebauer et al. (1987); (4) Guyon et al. (2006); 
(5) Glikman et al. (2006); (6) Sanders et al. (1989); (7) Elvis et al. (1994); (8) Neugebauer et al. (1979); (9) Balzano & Weedman (1981); (10) 
Rieke (1978); (11) McAlary et al. (1983); (12) Stein & Weedman (1976); (13) Surace et al. (2001); (14) Hayland & Allen (1982); (15) Gavazzi 
& Boselli (1996); (16) Veilleux et al. (2006); (17) Matsuoka et al. (2005); (18) Sitko et al. (1982); (19) Scoville et al. (2000); (20) Surace et 
al. (1999); (21) Imanishi et al. (2006); (22) Spinoglio et al. (1995); (23) Jarret et al. (2003); (24) Allen (1976) 



The high quality spectra made available by the IRS spec- 
trometer on-board the Spitzer Space Observatory (Houck et 
al. 2004) allows the detailed analysis of the MIR SED of a 
large number of AGNs. Here we fit the observed MIR spectra 
of 26 PG QSOs, already investigated by S08, using more real- 
istic three component models made of a clumpy torus, dusty 
NLR clouds and very hot dust clouds. In Sj2]we describe the 
observational data. In |J3] we detail our model and the fitting 
procedure and in |J4] we present the results of this procedure. 
Our main findings are summarized and discussed in |J5] 

2. SAMPLE SELECTION Spitzer OBSERVATIONS AND 
DATA REDUCTION 

Our sample consists of all AGNs in the QUEST Spitzer 
spectroscopy project (PID 3187, PI Veilleux). It is described 
in detail in Schweitzer et al. (2006; hereafter S06), Netzer 
et al. (2007) and S08. Most of the objects are Palomar- 
Green (PG) QSOs (Schmidt & Green 1983) and are taken 
from Guyon (2002) and Guyon et al. (2006). The luminosity 
range is Lsioo~ 10 44 5_46 ergs _1 where L5100 stands for XLx 
at rest wavelength 5 100A. Radio loudness and infrared excess 
are typical of these properties in low redshift PG QSOs. The 
QUEST sample, while representing the Guyon et al. (2006) 
sample, includes fewer sources and thus is not complete. This 
was explained in S06. Thus, several of the correlations dis- 
cussed below may differ from those in the entire sample in a 
way which is difficult to asses until a larger data set, with the 
same spectral coverage, is obtained. Here we aims at mod- 
eling the entire 2-35 fim SED hence we omit PG 1307+085 
(which appeared in the S06 sample) due to lack of observa- 
tions in the 5.2-8.7 /im range (Spitzer SL1 mode). Table Q] 
lists names, redshifts, and L5100 for all 26 QSOs in our sam- 



ple. 

The Spitzer observations and the data reduction of the 
QUEST sample are detailed in S06 and are summarized here 
for completion. IRS spectra for all objects were taken in both 
low resolution (5-14 fim) and high resolution (10-37 jj.m) 
modes. The standard slit widths of 3". 6 to 11". 1 include flux 
from the QSOs hosts and the vicinity of the AGNs. Data re- 
duction starts from the basic calibrated data (BCD) provided 
by the Spitzer pipeline. Specially developed IDL-based tools 
are then used for removing outlying values for individual pix- 
els and for sky subtraction. The SMART tool (Higdon et al. 
2004) is used for extraction of the final spectra. More details 
are given in S06. 

We have supplemented the Spitzer spectra with near in- 
frared (NIR) data obtained from the literature. The data are 
obtained from various sources, in particular Neugebauer et al. 
(1987) and the 2MASS extended and point source catalogs 
(Jarrett et al. 2000). For some of the sources there are several 
photometric measurements. For these, we average all mea- 
surements in each band and used the standard deviation as the 
flux error. For single observations, the flux error is estimated 
as 20%. These uncertainties must be underestimated since 
many of the NIR data were taken about 20 years before the 
Spitzer observations and all sources are variables. For a list of 
JHKL data references see table Q] 

We computed the bolometric luminosity of all sources us- 
ing their L5100 and a simple bolometric correction factor, BC, 
to convert it to Lbol- The best values of BC are discussed in 
Marconi etal. (2004), Netzer et al. (2007), and Netzer (2009). 
The approximation used here is taken from Netzer (2009) and 
is given by BC= 9 — logL44, where L44 =L5ioo/10 44 erg s _1 . 
For most objects the value of L5100 is obtained from the ob- 
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servations of Boroson & Green (1992). For PG 1244+026, 
PG 1001+054, and PG 0157+001 we use spectra from the 
SDSS (York et al. 2000) data release seven (DR7, Abaza- 
jian et al.2008) that were measured in the way described in 
Netzer and Trakhtenbrot (2007 ) . 

There are two uncertainties associated with the use of Lboi- 
The first is source variability, which is an important effect 
since the optical and MIR spectroscopy were separated by 
many years. We estimate this uncertainty to be a factor 
of ~ 1.5. The second uncertainty involves the approxima- 
tion used for BC. We estimate this uncertainty to be ~30%. 
Both uncertainties affect the derived model parameters such 
as torus covering factor and NLR distance. Because of this, 
we do not attach great importance to specific values in spe- 
cific sources obtained from our best acceptable models. On 
the other hand, the sample is large enough to enable a signif- 
icant analysis of its mean properties since the larger of these 
effects, due to source variability, is changing in a random way. 
The uncertainty due to the estimate of BC is smaller but more 
problematic since the expression we use can under-estimate 
or over-estimate Z^oi for all sources. This can introduce a 
systematic difference in torus and other component proper- 
ties, e.g. a smal ler c overing factors for all sources. We return 
to this point in ^5.11 

3. SPECTRAL MODELING 
3.1. Spectral Components 

The aim of the present work is to fit the 2-35 fim spectra 
of the sources in our sample. For this we combine models 
of three different physical components; a dusty torus, a dusty 
NLR and a very hot dust component. An additional starburst 
component representing emission from the host galaxy, is also 
present but is not part of the fit procedure. This component is 
subtracted from the observed spectra using the MIR spectrum 
of M82 in a way similar to the one used in S06. The reality 
of the subtraction can be judged from the remaining emission 
or absorption in the wavelengths corresponding to the strong 
PAH feature s (e. g. 6.2 or 7.7 fim). More details are given in 
S06 and in H3.2\ This section describes the spectral proper- 
ties of the three components and the following sections give 
detailed account of the fitting procedure. 

The first component represents a dusty torus surrounding 
the central energy source. AGN unification schemes suggest 
that the central engine is surrounded by dusty, optically thick, 
toroidal structure (e.g., Krolik & Begelman 1988; Antonucci 
1993). Earlier torus models assumed a smooth density dis- 
tribution (Pier & Krolik 1992; Rowan-Robinson 1995; van- 
Bemmel & Dullemond 2003). They all suffer from various 
incompletions and their agreement with MIR spectral obser- 
vations is generally poor. The more recent clumpy torus mod- 
els of N08 represent a significant improvement and are the 
basis for the present study. 

The fundamental difference between clumpy and smooth 
density distributions is that the former allows the radiation 
to propagate freely between different regions of the optically 
thick medium. The clumpy dust distribution results in the co- 
existence of clouds with a range of dust temperatures at the 
same distance from the central radiation source. They also al- 
low clouds at large distances to be exposed to the direct AGN 
continuum. In contrast, smooth density distribution models 
associate a certain dust temperature with a certain distance 
from the central source, thus limiting the range of acceptable 
SEDs. These differences allow the clumpy torus models to 



have a range of spectral properties not accessible in smooth 
density distribution torus models. N08 describe the formal- 
ism and the detailed radiative transfer used to calculate the 
MIR spectrum of clumpy dusty tori and our use of the torus 
models follows the same procedure. 

The first parameter of the clumpy torus model is the inner 
radius of the cloud distribution that is set to the dust sublima- 
tion radius Rj. This corresponds to a dust sublimation tem- 
perature r su b that depends on the grain properties and mixture. 
Here we adopt two different sublimation radii appropriate for 
two types of grains. The first is 
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where — Lboi/10 46 ergs -1 . This gives the innermost ra- 
dius where graphite dust can survive. It corresponds to a sub- 
limation temperature of about 1800K and is the one used in 
S08. N08 adopted a somewhat larger distance, which is more 
appropriate for silicate type grains with a sublimation temper- 
ature of 1500K. This is given by 
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(2) 



In the following we specify which of the two is used for each 
purpose. The clumpy torus model requires six additional pa- 
rameters: 

1 . The visual (5500A) dust optical depth of a single cloud, 
Ty (all clouds are assumed to have the same Ty). 

2. The mean number of clouds along a radial equatorial 
line, M). 

3. The ratio between the outer dimension of the torus and 
Rd, Y. 

4. The inclination angle of the torus with respect to the 
line of sight, i. 

5. The torus width parameter a, which is analogues to its 
opening angle. 

6. A parameter q that specifies the radial power-law dis- 
tribution of the clouds, i.e. N(r) r~ q , where N is the 
number of clouds. 

N08 also define an additional parameter that is set by the 
above free parameters. This is the probability P esc (fi) that 
light from the central source will escape the obscuring struc- 
ture at a given angle /3 without interacting with the clouds. 
Assuming that individual clouds are optically thick 



(3) 



where /3 = n/2 — i. By integrating P e sc over all angles and 
subtracting from 1, we obtain the probability of absorption by 
the torus. 



fi = l- 



k/2 



P eS c(j3)cos(j3) fl fj3. 



(4) 



This is also the fraction of obscured objects in a random sam- 
ple (e.g. the fraction of type-II AGNs out of the entire AGN 
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population), fj is equivalent to the "real" (geometrical) cov- 
ering factor of the torus, and does not depend on the inclina- 
tion angle since it reflects the global torus geometry. In other 
words, fi is the ratio between the total torus luminosity and 
Lbol- Due to the anisotropy of the torus radiation, fi differs 
from the apparent covering factor of the torus deduced from 
the ratio between its observed (angle dependent) luminosity 
and the bolometric luminosity of the central source. This ap- 
parent covering factor can be represented by 

1 riov, 

f(i) = — L(i,X)dX, (5) 

where L(i , A ) is the angle dependent monochromatic luminos- 
ity of the torus. For a more detailed description of all torus pa- 
rameters, as well as other model properties and assumptions, 
seeN08. 

The emission and absorption properties of the torus de- 
pend on the dust grains composition and other properties. All 
models consider here assume spherical dust grains with MRN 
size distribution (Mathis, Rumple & Nordsieck 1977). The 
dust has a standard galactic mix of 53% silicates and 47% 
graphite. The optical properties of graphite are taken from 
Draine (2003) and for silicates from Ossenkopf, Hennig & 
Mathis (1992) (OHM). The OHM dust mixture produces bet- 
ter agreement with observations of the 10 and 18 /im silicate 
features (Sirocky et al. 2008) and is the only one used in the 
present work. 

The second component of the model represents a collec- 
tion of dusty NLR clouds. The motivation for this component 
is explained in S08 where it was shown that such a compo- 
nent can contribute, significantly, to the MIR flux of luminous 
AGNs. The properties assumed here for these clouds are sim- 
ilar to the ones used in S08. We assume constant column 
density clouds with Nu = 10 21 ' 5 cm~ 2 . We further assume 
constant hydrogen density of 10 5 cm~ 3 , solar composition 
and galactic dust-to-gas ratio. The important physical param- 
eters for this component are the cloud-central source distance 
(which determines the dust temperature), the incident SED 
and the dust column density. The assumed gas density and 
column density only serve to define the emission from this 
component using convenient parameters such as the ioniza- 
tion parameter. The NLR component is assumed to be con- 
centrated in a thin spherical shell with a small covering frac- 
tion. These assumptions (that are reviewed later) mean that a 
single thin shell with the chosen properties is contributing to 
any 3-component model. 

The calculations of the emitted IR continuum of the clouds, 
as well as their emission line spectrum, are obtained with 
the photoionization code ION (Netzer 2006 and references 
therein). The assumed SED of the central continuum is the 
one described in S08 and the resulting NLR spectra are simi- 
lar to those shown in S08 Fig. 1. The clouds are optically thin 
to MIR radiation and hence no radiative transfer is required 
for the calculations of this emission. Obviously, this is not 
the case for the UV part of the spectrum where the opacity is 
much larger and where the transfer approximation we use is 
the one built into ION. 

The third component is a single BB representing hot dust 
emission. The need for such a component has been noticed 
in earlier works, e.g. in S08. This is very clear from the 
data, where a rise towards short wavelengths can be seen from 
the IRS spectra even without the supplemented NIR photom- 
etry. The simplified BB assumption is not entirely consistent 



with hot dust clouds illuminated from one side since for large 
dust optical depth, there must be a temperature gradient across 
such clouds. For small dust optical depth, the dust tempera- 
ture is more uniform but the cloud may be partly transparent to 
the incident radiation. For the purpose of the present study, we 
assume that all the incident optical-UV radiation is absorbed 
by this component and the dust te mpe ratur e is constant. We 
discuss this component further in £14.31 and fJ33]below. 

3.2. Fitting Method 

As shown in S08 and explained in detail in Netzer et al. 
(2007), the starburst contribution to the MIR spectrum can be 
significant, especially at long wavelengths. To remove this 
contribution, we follow the procedure of S08 who fitted and 
subtracted a "nominal" M82 spectrum from all spectra prior 
to the model fitting. We use the ISO-SWS mid-IR spectrum 
of M82 from Sturm et al. (2000) and the scaling factors of 
S08 that were obtained from the intensity of the strong PAH 
emission features that are clearly seen in some of our objects. 
Although this component is subtracted prior to the model fit- 
ting, its normalization introduces another degree of freedom 
to the procedure. The remaining spectra are assumed to rep- 
resent the intrinsic AGN continuum. 

The fitting procedure minimizes the modified x 2 , 

2 1 /vL v ,obs(A) - VL v ,model(A) \ 2 

X =^{— a x J ' (6) 

where A^of is the number of degrees of freedom, vL v , bs(A) 
is the starburst subtracted luminosity deduced from the ob- 
served monochromatic flux density at rest wavelength A, and 
model (A ) is the combined emission due to all three com- 
ponents, 

^v, model 

(A) = L v ,toras(A) +flNLRiv,NLR(A) + aBB^V.BB (A ) , 

(7) 

where «nlr and abb are fitting parameters. The bolometric 
luminosity of the source uniquely determine the total emer- 
gent flux of the torus component through the radiative transfer 
calculation. Since an independent measure of the bolometric 
luminosity is available for these type-I sources, the normal- 
ization of the torus component is not a free parameter of the 
fitting procedure (i.e. the first term in Eq.|7]is set by the torus 
parameters but the total energy emitted by the torus is fixed. 

To reduce the uncertainty, we first bin the observed spectra 
into ^100 bins of equal energy widths and take in each the 
mean observed flux. The torus, NLR and BB models are then 
interpolated to the same energy grid. The binning of the data 
practically removes all emission lines and smooth the spec- 
trum, allowing to focus on the broad band continuum shape. 
Some remaining peaks may result in larger local % 2 values but 
since none of the models contains lines, the effect would be 
the same for all fits. 

According to S08, most of the errors in the IRS spectra 
are systematic and are approximately 10% of the flux. We 
add this systematic error to the standard deviation of the data 
points in each energy bin to get the value of Ox in Eq. [6] 
This is taken to represent the error on the observed binned 
flux. The use of this a affects the derived % 2 sucn that its 
absolute value is different from the one normally used in sta- 
tistical analysis. However, the relative value of % 2 can still 
be used to discriminate between different models and will be 
treated as such in the remaining of the paper. 
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Rest wavelength (Jim) Rest wavelength dim) 

FIG. 1 . — Best torus-only fit to the spectrum of PG 1004+130 (left; best fit of this type) and PG 1440+356 (right; typical torus-only fit). The observed binned 
data are shown in black, asterisks represent K, L photometry, and the best fit model is shown in red (top panels). The quality of the fit is demonstrated by the ratio 
between the data and the model (middle panels) and the Rvalue in each wavelength bin (bottom panels). Note the clear flux deficiencies at short wavelengths in 
both fits and the additional flux deficiencies at longer wavelengths for PG 1440+356 



The fitting algorithm computes % 2 values for all possible 
combinations of torus, NLR and BB models. There are six 
free parameters in the tor us m odel that describe its geometri- 
cal properties (see section [3j]i. The radial extent of the torus, 
Y, ranges from 5 to 60 R c i. The total number of clouds along 
radial equatorial lines ranges from 1 to 10. The optical depth 
of individual cloud, Ty, ranges from 5 to 100. The torus width 
parameter, <j, ranges from 15° to 75°. The index of the power 
law radial distribution of clouds, q, only accepts one of the 
three integer values, 0, 1 or 2. 

Finally, the torus inclination angle, i, ranges from 0° to 90°. 
In the NLR model there are two free parameters, the cloud 
distance and the normalization factor «nlr- The distance is 
changed in steps of 0.075 dex between 1 pc to 850 pc for a 
source with Lbol = 10 45 ergs _1 . For the hot BB, there are also 
two free parameters, temperature and a normalization factor 
«bb- The range in temperatures is 800-1800 K. 

4. RESULTS 

4.1. Torus only models 

We first examine the possibility that single torus models can 
explain the observed spectrum. This is important in order to 
evaluate the general properties of clumpy tori models and has 
never been attempted on large Spitzer -IRS samples. In gen- 
eral, allowing for a torus-only model results in large % 2 and 
poor fit to almost all spectra. On average, % 2 values of such 
fits are about an order of magnitude larger than those obtained 
by using all three components. None of the torus models 
could fit well the full observed spectrum of any of the objects. 
The best torus-only models are those fitted to the spectra of 
PG 1626+554 ,PG 1004+130 and PG 0838+770. These ob- 
jects are fitted reasonably well around the silicate emission 
features and at long wavelengths but fail at A < 8jtim where 
they show clear flux deficiencies (see left panel of Fig.Q]). Sin- 
gle torus fits to all other objects are significantly worse. This 
is illustrated in the right panel of Fig. Q] that shows the best 
torus-only fit to the spectrum of PG 1440+356. Not only that 
the model poorly fits the silicate features region, it all shows 
insufficient flux both at shorter and longer wavelengths. 



4.2. Combined torus-NLR models 

Next we examine the quality of two component models. In 
general, the fits improved compared with torus-only models 
especially over the wavelength range A >6 jim. An example 
is shown in Fig. [2] where we reproduce the best two compo- 
nent fit to the spectrum of PG 141 1+442. In this case, the fit 
is within 10% of the observed flux for all A >5 jxm. Clear 
deviations are still found at the short wavelength part of the 
spectrum. We were able to obtain satisfactory fits over the 
^6-35 jJ.m range for 17 out of the 26 objects using two com- 
ponent models. The x 2 values for these fits are, on average, a 
factor 6 smaller than the single torus fit values. These results 
are comparable to those shown in S08 where a combination 
of several BBs and NLR components was fitted to the 5-35 
jj,m spectrum of all sources. Only in PG 1626+554, which is 
not part of the above group of 17 objects, the torus-only fit is 
superior to the two-component fit (i.e. the NLR component is 
not required). For the remaining 8 objects, the % 2 values of 
the two-component fits are only about a factor 2 smaller than 
the values obtained for torus-only fits. In these 8 fits, none of 
the main spectral features are well fitted. 

In all cases fitted with by two component models, a signifi- 
cant fraction of the flux shorter of ^6 fim is still unaccounted 
for. The mean missing fractional (i.e. the difference in flux 
between model and observations over the 2-35 jJ.m range) is 
about 40%. 

4.3. Three component models 

The missing flux in the short wavelength part of the spec- 
trum in most of our sources is the main motivation for our 
third, hot dust component. This was not discussed in S08 who 
fitted only the A > 5/im part of the spectrum. 

We have investigated two possibilities: The first is emis- 
sion due to hot dust with the same grain composition used 
for the other components. This can represent additional dusty 
clouds in the innermost part of the torus that are not part of 
the clumpy structure of N08. Including such a component, 
obtained from the single cloud models of N08, resulted in 
poorer agreement between model and observations. The rea- 
son is that such a component produces strong silicate emission 
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Rest wavelength (pm) 

FIG. 2. — Best fit to the spectrum of PG 141 1+442 (red line) using a combi- 
nation of torus (dashed line) and NLR (thick black line) components. Middle 
and bottom panels are as in Fig[T] 



TABLE 2 

X 2 FIT RESULTS 



Object 


Torus-only fit 


2-component fit 


3-component fit 


PG 2349-014 


6.80 


2.34 


0.14 


PG 2251+113 


8.99 


3.95 


0.25 


PG 2214+139 


2.63 


1.99 


0.27 


PG 1700+518 


2.03 


0.77 


0.24 


PG 1626+554 


1.10 


1.12 


0.50 


PG 1617+175 


3.24 


1.14 


0.10 


PG 1613+658 


2.62 


1.00 


0.12 


PG 1448+273 


1.56 


1.24 


0.31 


PG 1440+356 


3.36 


1.34 


0.15 


PG 1435-067 


3.59 


0.52 


0.31 


PG 1426+015 


2.15 


0.93 


0.11 


PG 1411+442 


2.14 


0.93 


0.10 


PG 1309+355 


3.21 


1.52 


0.88 


PG 1302-102 


3.08 


0.70 


0.10 


PG 1244+026 


1.03 


0.49 


0.34 


PG 1229+204 


2.45 


0.94 


0.18 


PG 1 126-041 


4.86 


1.51 


0.24 


PG 1116+215 


5.87 


2.74 


0.24 


PG 1004+130 


1.90 


0.58 


0.23 


PG 1001+054 


7.91 


2.23 


0.28 


PG 0953+414 


7.01 


4.15 


0.45 


PG 0838+770 


1.61 


0.79 


0.26 


PG 0026+129 


7.62 


2.07 


0.31 


PG 0157+001 


12.40 


1.80 


0.81 


PG 0050+124 


8.26 


1.39 


0.20 


B2 2201+31A 


4.61 


1.19 


0.09 



features, especially the broad 10 fim feature, that completely 
change the SED. The apparent strength of the 10 jim and the 
18 jJ.m features in pure emission models represents the aver- 
age temperature of the dust. The introduction of very hot sili- 
cate dust gives the impression of a very strong 10 jJ.m feature 
compared with a weak 18 /im feature, which clearly contra- 
dicts the observations. Since the hot dust luminosity is very 
significant, the variations in the intensities of the other com- 
ponents cannot compensate for this increase. Thus, using hot 
silicate dust we could not successfully fit the observations. 

The second possibility is a hot BB. The inclusion of the this 
component improves, significantly, the quality of most fits. It 
allows for a better fit at both short and long wavelengths since 
the torus and NLR components are no longer needed to ac- 
count for the short wavelength part of the spectrum. For the 
same reason, the models come closer to the observations at 
10 jxm and 18 jJ.m. The % 2 values obtained with the three- 
component fits are, on average, an order of magnitude smaller 



than the torus-only fits and about a factor of 3 smaller com- 
pared with the two-component fits. % 2 values for all objects 
and all types of fits are listed in table [2] (note again that these 
values do not represent the formal % 2 values). Given all three 
components, we obtain fits that are within 20% of the obser- 
vations over the entire wavelength range. In the wavelength 
range of the silicate emission features (^8 — 20 fim), 19 of 
our best fit models deviates by no more than 5% from the 
data. In the remaining 7 sources, the deviation over this wave- 
length range is at most 10%. A more realistic model for the 
short wavelength com ponent is a collection of hot graphite- 
only clouds (see i!5.3l ). Such clouds have no silicate features 
in their spectra. Model fitting including such a component is 
deferred to a future publication. 

Fig. [3] shows the best fit three-component models for two 
representative cases, B2 2201+31A and PG 1617+175. The 
top panel of each diagram shows the best fit model (in red) and 
the observed spectrum (in black). We also show individual 
components: torus (dashed line), NLR (thick black line) and 
hot BB (dash-dotted line). Asterisks represent the K and L 
photometry. In the bottom and middle panels of each diagram 
we show the quality of the fit in each wavelength bin. 

Having obtained satisfactory fits for all sources, we have 
calculated the median contribution, over the entire sample, of 
each of the components. This is done at every wavelength and 
is shown in Fig. [4] in sucn a way that the sum of all compo- 
nents at each wavelength bin is 1 . The diagram shows that the 
BB component dominates the spectrum below ^4 jim while 
the torus dominates between 5-25 fim. The NLR component 
has a non-negligible contribution (~ 40%) above ^16 jim, 
including the 18 jlim emission feature. 

4.4. Fit quality 

Each of the fits is assigned a quality flag set by four criteria, 
the calculated % 2 over the entire 2-35 jim range and partial 
X 2 values for three shorter segments of the spectrum. Re- 
garding the first, since the errors are not normally distributed, 
the best % 2 has a different meaning from the commonly used 
X 2 and setting a standard confidence level is not applicable. 
Given this, we decided to consider only fits with % 2 values 
not exceeding the minimum % 2 by more than a factor of 1.25. 
This conservative narrow range of acceptable % 2 was chosen 
to ensure that it contains, beside the best fit, only acceptable 
fits. Indeed all values within this range seem very good and in 
most cases, there are also a few acceptable fits outside of this 
range. However, the distribution of the model parameters for 
the acceptable fits is not affected much by these outliers. 

The three additional partial J 2 criteria are aimed to identify 
those fits where a certain part of the spectrum is not fitted well 
even though the global % 2 is within the above mention range. 
For this, we examine the goodness of the fit in three different 
parts of the spectrum by calculating % 2 values for the ranges 
~ 2 — 7 jj.m, 7—20 ftm and ~ 20 — 35 jJ.m. We used a similar 
approach (in this case a factor of 1.5 relative to the best x 2 ) 
and reject fits that are outside of this range. 

Having fixed the % 2 limits, we selected the fits that pass all 
four criteria (one global and three partial % 2 limits). There are 
thousands of possible torus models and a much larger number 
of possible combinations hence, the number of accepted fits 
is large, typically 10-50 per sources. All these are inspected 
visually to verify the success of the automatic selection pro- 
cess. The groups of acceptable models, for all sources, is the 
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FIG. 3. — Three component fits to the spectra of B2 2201+31A (left) and PG 1617+175 (right). The addition of a BB component (dash-dotted line) clearly 
improve the quality of the fit for these objects and the rest of our sample. Middle and bottom panels are as in Fig[T] 



basis for the following analysis. 

Two objects in our sample, PG 0157+001 and 
PG 0050+124, have total covering factors larger than 
unity. This can perhaps be explained by source variability. 
To check this for the entire sample we compare the Boroson 
& Green (1992) data with more recent observations from the 
SDSS. New spectra are available for 8 objects in our sample, 
5 of which exhibit an increase in luminosity by a factor ~ 2, 
two show a decrease of about the same factor and one shows 
no significant change. Of the two objects with a covering 
factor larger than unity, only PG 0157+001 has recent SDSS 
spectra. A comparison with the older observations show no 
significant change in luminosity between the two epochs. 
Under-estimation of the star formation contribution to the 
mid-infrared range of the spectrum may also result in an over 
estimation of the covering factor. Obviously, we could not fit 
the spectra of these two objects with adequate quality, using 
the measured value of their bolometric luminosity (see the 
discussion rega rding the implications of Lbol related uncer- 
tainties in ^5! lb . To proceed, we assume that flux variation 
is the cause for the apparent discrepancy and artificially 
lower the bolometric luminosity of these objects so that their 
covering factors become 1 . This allows us to fit their spectra 
in a way similar to the other sources and get adequate fits and 
reasonable model parameters. However, we do not include 
them in the overall statistics and all parameter distributions 
and mean values presented here do not include these objects. 

5. DISCUSSION 

Having found the best three-component models for all indi- 
vidual spectra, we now discuss the properties of these compo- 
nents as well as their distribution in the sample. 

5.1. Torus properties 

We examine the distribution of the different torus parame- 
ters for each of the objects in our sample. Some of the pa- 
rameters have a narrow distribution around a mean value and 
in some cases only a single acceptable value. Other parame- 
ters exhibit a broader, more uniform distribution. In Fig.[5]we 
present a typical example of the acceptable torus parameters 
for one source, PG 2349-014. As seen in this case, C, i, and q 
have very narrow distributions for all acceptable torus models. 
No ,Y, and Ty have broader distributions around their mean 
values. The values obtained for individual objects from their 
best fit models are listed in table|3]The parameter distributions 




10 12 14 16 18 20 22 24 
Rest wavelength (urn) 

FIG. 4. — Relative contributions of the BB component (dash-dotted line), 
torus (dashed line), and NLR (solid line) to the fitted models (median values 
at each wavelength for the entire sample). 
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FIG. 5. — Typical example of the distribution of acceptable torus parame- 
ters. The objects is PG 2349-014 and the arrows denote mean values. Note 
that O", q and i have a very narrow distribution while No, Y and Ty have a 
broader one. 



of all sources were combined by giving each value within the 
acceptable range in a certain source its relative weight in the 
distribution. The results of this combination are general pa- 
rameter distributions that are shown in Fig[6] Inspection of the 
various parts of this figure lead to the following conclusions: 
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TABLE 3 

Torus parameters for best fit models 



Object 


a 

(deg) 


Y 


N 


1 


T V 


(deg) 


FLr 2349-014 


29 


35 


7 


1 


80 


48 


FLr 2251 + 11 j 


15 


23 


S 


1 


49 


67 


PU 2214+139 


25 


25 


4 


1 


7 1 


40 


PCj 1 700+5 Is 


43 


34 


3 


2 


87 


43 


PU 1626+554 


28 


6 


6 


1 


17 


31 


FLj 161 /+175 


46 


44 


2 


2 


93 


20 


FLr 1o1j+65s 


48 


55 


3 


1 


80 


27 


FCj 144o+27j 


19 


1 1 


X 


1 


59 


53 


PU 1440+356 


40 


12 


4 





41 


26 


P(j 1435-06/ 


57 


34 


1 


2 


7 1 


38 


PCj 1426+015 


35 


33 


5 


1 


75 


30 


PCj 141 1+442 


55 


40 


2 


2 


70 


14 


PCj 1309+355 


22 


24 


6 


1 


67 


37 


tic i i m 
PCj 1 302-102 


37 


30 


6 





44 


32 


PCj 1244+026 


33 


40 


6 


1 


68 


31 


PG 1229+204 


34 


33 


5 


2 


91 


28 


PG 1 126-041 


34 


36 


6 


1 


61 


28 


PG 1116+215 


34 


31 


3 


1 


64 


39 


PG 1004+130 


31 


39 


4 





38 


35 


PG 1001+054 


22 


29 


7 


1 


72 


38 


PG 0953+414 


18 


20 


7 


1 


66 


35 


PG 0838+770 


49 


28 


4 





62 


23 


PG 0026+129 


27 


35 


3 





79 


43 


PG 0157+001 


48 


52 


5 





71 


16 


PG 0050+124 


51 


35 


5 





20 


8 


B2 2201+31A 


16 


33 


8 


2 


35 


29 



1. The mean torus width parameter is a — 34°. None 
of the fits requires the largest allowed a of 75°. This 
value determines, more than any of the other parame- 
ters, the geometrical covering factor (J2) of the torus 
and is consistent with the expected type-I:type-II ratio 
(see below). 

2. The mean radial extent of the torus, Y, has a broad dis- 
tribution with a mean value of 3 1 . The range in Y im- 
plies torus outer sizes, which range from ^1 to 35 pc 
using R c i c and ^3 to 90 pc using RdSi- 

3. The average number of clouds along an equatorial ray, 
No, has a broad distribution with a mean of 5 clouds. 
The distribution is practically uniform for No>3. 

4. The parameter that specifies the radial power-law dis- 
tribution of the clouds, q, has a mean value of 1 . This 
parameter is related to the anisotropy of the torus ra- 
diation. As q decreases the torus radiation becomes 
less isotropic. The SED shape is affected less by q and 
hence we do not attach great significance to its specific 
value. 

5. The mean optical depth of a cloud is Ty=58. This pa- 
rameter, again, has a broad distribution covering the 
range Ty > 20. Since the requirement for large MIR op- 
tical depth is Ty-^ 10, it is not surprising that the torus 
models are not very sensitive to the exact value of Ty 
beyond this value. 

6. The torus inclination angle shows a broad distribution 
for all i < 60° with a mean of 33°. This again is consis- 
tent with the assumption that the direct line-of-sight to 
the center of type-I AGNs is almost completely free of 
obscuring material. 
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FIG. 6. — Torus parameter distributions for the sample excluding 
PG 0157+001 and PG 0050+124 (see text; arrows denote mean values). 



We also checked for various correlations between those pa- 
rameters and found no such correlation. This suggest that 
torus models with different parameters can result in similar 
SEDs. 

As explained, SED related uncertainties, in particular the 
bolometric luminosity of the sources, make it difficult to de- 
termine the exact torus properties in individual objects. There 
is a relatively large range in some of the parameters (e.g. No 
and i) and small changes in Lj,oi can result in a significant 
change in those parameters. For example, lower bolometric 
luminosity, due to changes in L5100 between the optical and 
the Spitzer observations, requires lower L Vj torus and higher in- 
clination angle. In principle this may introduce non-negligible 
uncertainties on ; (and other parameters) in individual sources. 
However, the average sample properties are less likely to de- 
pend on such uncertainties and the discussion below focus on 
this aspect of the work. Below we examine three of the param- 
eters, i, Nq, and a using Eq.[3] We set two of the parameters 
to their mean values and investigate the dependence of P e sc on 
variations of the remaining free parameter. This is equivalent 
to the investigation of a generic torus that represents the entire 
sample of 26 AGNs. 

First we consider the torus inclination angle i as the free 
parameter. Using Eqf3] together with the mean values found 
here, the probability that light from the central region will 
escape such a structure, and the source will be classified as 
type-I AGN, is about 75%. Setting No and a to their mean 
values (5 and 34°, respectively) and changing the inclination 
angle, we find that for i = 50° P esc ~ 30% and for i = 70° 
Pesc < 3%. We therefore expect that for type-I objects, the 
inclination angle should not drop below ~ 60°. As can be 
seen in Fig|6] the distribution of the inclination angle is indeed 
consistent with random inclination angles in the range — 60° 
with no preference to a specific angle. 

The second important parameter is the mean number of 
clouds, jVo- Setting a and i to their mean values and changing 
M), we find that P e sc remains large even for very high values 
of No- Thus, the exact value of No is not very important pro- 
vided the inclination angle is not too large. This is the reason 
why No has the broad distribution seen in Fig|6] 

Next we tested the acceptable range of the torus width pa- 
rameter, c, by fixing all other parameters to their mean values. 
We find that P e sc is sensitive to such changes and falls rapidly 
for a ^45°. This is the reason that the a distribution falls 
rapidly for a > 45°. Torus with high value of a would have 
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TABLE 4 

rMDDCI ATTHW D \7AT TTT7C 

LUKKnLAl 1UJN r VALUED 












Pearson 




Spearman 




Parameter 


log(L bol ) 


log(M Bli ) L/Lem 


/og(Lboi) 


log(M BH ) 


L/LEdd 


a 
Y 


0.031(-) 




0.028(-) 






No 
9 
Tv 
i 


0.014(+) 




0.0047(+) 




... 


^esc 
^NLR 

BB temperature 
Real total CF 
Apparent total CF 

h 

m 

NLR CF 
BB CF 
/ 2 +BB CF 

A^torus 


0.00018(+) 

0.0033(-) 
0.001 1(-) 
0.026(-) 
0.005 1(-) 

0.013(-) 
0.0046(-) 
0.013(+) 


Z 0.0054(+) 
0.01(+) 

0.024(-) 
0.0063(-) 

0.026(-) 

0.048(-) 
0.019(+) 


0.0004(+) 

0.0056(-) 
0.0013(-) 
0.029(-) 
0.0026(-) 

0.04(-) 
0.0045(-) 
1.7e-06(+) 


0.045(-) 
0.014(-) 

0.016(-) 
0.0003(+) 





M torm /M m 

Morus/^bol 



0.0083(+) 



NOTE. — Plus and minus signs indicate positive and negative correlations, respectively. Results for p > 0.05 are not listed. 



low P esc and thus a low probability of being classified as type-I 
AGN even if the inclination angle is zero (a face-on torus). 

The combined effect of the above parameters is, perhaps, a 
more meaningful test. In particular, the combination of a and 
Nq determine the geometrical covering factor of the torus, fa, 
which is combined, later, with the equivalent prop erty of the 
other components. This issue is discussed in ij5.4| 

The above procedure is not applicable to the visual optical 
depth parameter. This parameter does not affect f esc or fa but 
influence the SED shape, in particular the apparent strength of 
the silicate emission features. The mean visual optical depth 
for the entire sample is Ty=58 and the lowest value is 20. 
Thus, realistic clumpy torus models require large IR optical 
depths to explain the observed spectrum. The exact value of 
Ty is not very important once the clouds are thick enough. 

The torus size found here is consistent with upper limits 
set by several high resolution observations of nearby AGNs 
(e.g, Soifer et al. 2003; Radomski et al. 2003; Jaffe et al. 
2004). This is in contrast to previous theoretical works involv- 
ing smooth density distribution torus models. In these models, 
the torus is required to be much larger (up to few hundreds of 
pc) in order to posses the large amounts of cool dust neces- 
sary for producing the observed MIR emission. This require- 
ment arises because in smooth density distribution models, 
the dust temperature is uniquely determined by the distance 
from the center. In contrast, the clumpy torus model enables 
emission from cool dust clouds that are much closer to the 
center thus able of producing the required MIR emission with 
much smaller dimensions. 

Perhaps the more interesting results are the correlations of 
the torus parameters with the physical properties of the AGNs, 
in particular BH mass (MbhX bolometric luminosity and nor- 
malized accretion rate, L/L^dd- To study these correlations 
we have estimated these quantities for all sources using the 
procedure described in Netzer & Trakhtenbrot (2007). In this 
procedure (the "virial" mass determination, see Eq. 1 in Net- 
zer & Trakhtenbrot 2007), L 5100 and FWHM(Hj3) are com- 
bined to obtain Mbh, and L/LEdd is obtained by using the as- 



sumed bolometric correction factor BC. Table|4]lists p-values 
for Pearson and Spearman-rank correlation tests for the dif- 
ferent torus parameters against those properties. These val- 
ues represent the probability that the observed correlation oc- 
curred by chance. We regard as significant all correlations 
where p < 0.05 and indicate with plus and minus signs posi- 
tive and negative correlations, respectively. Out of the above 
geometrical parameters, only a and i show a significant cor- 
relation with Lbol- This is consistent with a 'r eced ing torus' 
assumption (Lawrence 1991) and discussed in ; 



5.2. NLR properties 

The important parameters of the NLR component are its 
distance from the center, which determines the dust tempera- 
ture, and its covering factor, which determines the fraction of 
MIR flux emitted by the clouds (anlr m Eq|7]). 

The photoionization models used to calculate the NLR 
emission are all of the same constant density («h = 10 5 cm -3 ) 
and column density (10 215 cm -2 ). As explained earlier, the 
clouds are optically thin in the MIR range and their distance 
from the central source uniquely determine the NLR dust tem- 
perature. The assumed gas density plays no important role 
and serves only as a convenient way to determine the pho- 
toionization model parameters (in the actual calculations, the 
model parameters are density and ionization parameter). The 
dust column density is determined by the assumed gas com- 
position and depletion. The present calculations use solar 
composition with abundances smaller than assumed in sev- 
eral recent NLR models such as the constant pressure models 
of Groves et al. (2006). This again is of little importance since 
the dust column density is the only important parameter, un- 
less an unusual grain mixture is used. 

Table [5] lists NLR distances for the 26 best models. Radii 
are given in pc as well as in units of the sublimation distances, 
Rd C an d Rd Si- F° r the graphite sublimation radius, the NLR 
distances span the range 256 — 1755/?^,c with a mean value of 
707 Rdc- This translates to 101 — 694Rdsi with a mean value 
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FIG. 7. — NLR distance distribution in units of the graphite sublimation 
radius. 



in' - 




FIG. 8. — NLR cloud distances against Lb„i. Circles represent the NLR 
distances in pc, and the solid line is a linear fit to the data with a slope of 
0.47. The graphite dust sublimation radius is shown for comparison as a 
dashed line. 

of 280R ( / si- The distribution of ^nlrIH units of the graphite 
sublimation radius is shown in figure [7] 

The mean distance found here is about a factor of 4 larger 
than the values found by S08. The reason for the larger dis- 
tances is mostly due to the introduction of real torus mod- 
els compared with the combination of BBs used in S08. The 
clumpy torus models contribute a certain amount to the ob- 
served silicate emission and hence the NLR contribution, 
which was the only one with such features considered in S08, 
is reduced. This can also be seen when comparing the NLR 
covering factor found here and the one found in S08 (see Ta- 
ble |6). The mean covering factor for the NLR component in 
S08 is a factor 2 larger than the one found here which is ~7%. 
We consider the prese nt values more appropriate to the objects 
in question (see ^5.4b 

Figure[5]shows NLR cloud distances against the bolometric 
luminosity of the objects. Using the values in Table[3]we find 
the following scaling relationship, 

fl NLR = 295 x1°™ 13 pc- (8) 

The slope is consistent with 0.5, which is the one predicted for 
simple single-cloud photoionization models and a constant U 
gas. 

The dusty NLR dimensions can be compared with earlier 
works on NLR size such as Bennert et al. (2002), Schmitt 
et al. (2003), Netzer et al. (2004) and Baskin and Laor 



TABLE 5 

NLR PROPERTIES FOR BEST FIT MODELS 



Object 


log(f/f 


^NLR 
(PC) 




^NLR/^rf,S( 


PG 2349-014 


-2.6 


201 


608 


240 


PG 2251+113 


-2.2 


323 


400 


158 


PG 2214+139 


-2.0 


58 


271 


107 


PG 1700+518 


-2.9 


780 


913 


361 


PG 1626+554 


-1.6 


161 


726 


287 


PG 1617+175 


-2.8 


157 


832 


329 


PG 1613+658 


-1.6 


152 


518 


205 


PG 1448+273 


-2.2 


74 


392 


155 


PG 1440+356 


-3.5 


308 


1755 


694 


PG 1435-067 


-2.7 


155 


737 


291 


PG 1426+015 


-1.8 


57 


256 


101 


PG 1411+442 


-2.7 


135 


701 


277 


PG 1309+355 


-2.1 


161 


363 


144 


PG 1302-102 


-2.2 


210 


431 


170 


PG 1244+026 


-2.8 


79 


867 


343 


PG 1229+204 


-2.5 


69 


558 


221 


PG 1126-041 


-3.0 


109 


946 


374 


PG 1116+215 


-2.6 


270 


579 


229 


PG 1004+130 


-2.3 


302 


470 


186 


PG 1001+054 


-2.5 


151 


603 


239 


PG 0953+414 


-2.0 


145 


318 


126 


PG 0838+770 


-2.0 


57 


345 


136 


PG 0026+129 


-2.4 


120 


428 


169 


PG 0157+001 


-3.4 


597 


1520 


601 


PG 0050+124 


-3.6 


417 


1728 


683 


B2 2201+31A 


-3.1 


1234 


1121 


443 


Mean 




249 


707 


280 



a Assuming n = 10 5 cm 



(2005). These papers use two different methods to estimate 
the NLR size. The Bennert et al. and Schmitt et al. works are 
based on detailed HST imaging of AGN samples that are of 
the same redshift range as the one considered here. Our dust 
temperature-based dimensions are smaller by a factor of ~ 10 
compared with these imaging-based estimates. This is not 
surprising since our estimates are biased toward those NLR 
clouds that have the largest covering factor for the hottest 
dust clouds. The discrepancy suggests that real NLRs, i.e. 
those with a range of dust temperatures contribute more to the 
longer wavelength part of the spectrum. Furthermore, Netzer 
et al. (2004) suggested that the Bennert et al. (2002) estimates 
are considerably larger than "real" NLR sizes for two reasons. 

1/2 

The first is that the 7?nlr 06 ^; on dependence must break down 
at some intermediate luminosity since extrapolating this rela- 
tion to high luminosity objects would result in NLR dimen- 
sions that exceed the host galaxy dimension. The second rea- 
son is that Bennert et al. (2002) measured image sizes that 
encompasses more than 98% of the detectable [Olll] A5007 
emission. A 90% or 95% intensity limits can lead to much 
smaller dimensions. 

Baskin & Laor (2005) used a different method to estimate 
the distances to 40 NLRs in the Boroson & Green (1992) sam- 
ple. They fitted the observed intensities of the [Olll] A 5007, 
[O III] A 4363 and Hj3 narrow emission lines using two differ- 
ent models. The first, assumed a single uniform, [O III] A5007 
emitting region. This resulted in ^nlr ~ 130 x L°g 45 pc 3 and 
a very large mean gas density of 10 5,85 cm~ 3 . These dimen- 
sions are smaller than those found by Bennert et al. (2002). 
The second model assumed two emitting regions: a low, con- 

3 In this expression and the following ones, we converted the continuum 
luminosity given by Baskin & Laor (2005) into bolometric luminosity. 
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stant gas density (« = 10 3 cm 3 ) region where most of the 
[O III] A5007 originates and a high, constant gas density (n = 
10 7 cm~ 3 ) region where most of the [Olll] A4363 originates. 
For the low density region they found /?nlr ~ 1850 x L^ 34 pc 
and for the high density region /?nlr ~ 4 x L^ 5 pc. The 
[Olll] A5007 dimensions are comparable to those of Ben- 
nert et al. (2002). It is clear that single-zone and even two- 
zone NLR models are highly simplified. Moreover, dimen- 
sion deduced from a narrow line intensity is uncertain since 
the gas density is changing across the NLR and in general 
/?nlr x n~ ' 5 . Thus, it is not surprising that the NLR size 
obtained from such estimates spans a very large range. 

The M82 starburst template subtracted from the spectra 
could, in principle, have an effect on the calculated NLR prop- 
erties, e.g. a different template could have larger contribution 
towards shorter wavelengths (^30 fim). Consequently, the 
NLR component would have less weight and different spec- 
tral shape that would result in larger NLR distances. We re- 
gard this possibility as an additional uncertainty on the deter- 
mination of the NLR distance. 

In conclusion, our hot dust measure of the NLR size is 
based on a method different than the direct imaging method 
and the photoionization modelling method. While further 
discussion of those discrepancies is beyond the scope of the 
present paper, the main conclusion is that simple, dusty NLR 
clouds (or single NLR shells) at the distances found above, in 
combination with the two other model components, can ade- 
quately fit the 2-35 jum spectrum of the QUEST AGNs. 

5.3. Hot BB properties 

The 2 — 4jim wavelength range is dominated by the hot BB 
component of our model (see Fig. 0). On average, the lumi- 
nosity of this component is about 40% of the total 2 — 35/im 
luminosity and is comparable to the luminosity radiated by 
the clumpy torus. The mean BB temperature is 1400 K and 
the distribution of temperatures is shown in Fig. [9] While 
such a large contribution is evident in all of the sources, the 
physical origin of this component is not yet clear. Clumpy 
tori of the type considered here cannot produce more than one 
luminosity bump. Since most of the torus radiation peaks at 
5-20 jim, the short wavelength excess cannot result from this 
structure. The conclusion is that the assumed BB must be a 
separate component. 

We have considered several possible explanations for the 
hot BB component. The first is a contribution from old stellar 
population near the galactic center. To examine this, we use 
Mgn to estimate the host galaxy luminosity using the scaling 
relations of Lauer et al. (2007; see their Eq. 3). To convert 
visual magnitudes to K-band magnitudes we use (V-K)=3.2 
typical of giant elliptical galaxies (e.g. Grasdalen 1980). This 
K-band luminosity was compared to the luminosity of the BB 
found here. In all cases, the BB component is about a factor 
10 more luminous than the derived K-band luminosity of the 
host galaxy. Thus, such a stellar population would not have 
sufficient luminosity to account for the observed 2-7 jj,m flux. 

The space distribution of such an old stellar population is 
another limitation. Broad band monitoring of several nearby, 
intermediate luminosity AGN clearly show K-band variations 
on time scales of weeks to months (e.g., Suganuma et al. 
2006). This has been interpreted as changes in the tempera- 
ture of the innermost dust, just beyond the sublimation radius. 
Given the great similarity between the 2-35 jim spectrum of 
the present sample and the spectra of those variable AGNs, 
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Black body temperature [K] 
FIG. 9. — Blackbody temperature distribution for the hot dust component. 

we conclude that the origin of the K and the L-band emission, 
in the present sample and the above AGN sample, is the same. 
Such dimensions are clearly too small for the assumed stellar 
population. Needless to say, any K or L-band variations are 
inconsistent with stellar population properties. 

Another strong source of radiation is starburst activity in 
the host galaxy. The IR emission from starburst, however, is 
expected to peak at much larger wavelengths. S06 found a 
strong correlation between the strength of the PAH feature at 
7.7 fim and the continuum luminosity at 60 jJ.m. The same 
trend is found in starburst dominated ultraluminous infrared 
galaxies (ULIRGs; see Fig. 4 in S08). Thus, luminosity pro- 
duced by starburst in the QUEST sample dominates the far-IR 
part of the spectrum and cannot account for the hot dust com- 
ponent. 

Another possibility that has already been mentioned is that 
the new component represents emission by very hot dust, 
wh ich i s not a part of the clumpy torus structure. As explained 
in H4.3\ we could not fit the observations assuming single hot 
dust clouds of the same grain composition as the torus. The 
reason is that such clouds produce extremely strong 10 /im 
emission. This is not the case for pure graphite dust that has 
no prominent features at MIR wavelengths. Graphite dust 
grains can survive at the BB temperatures found here and pro- 
vide a suitable explanation for the BB component. The loca- 
tion of this dust can be inside the silicate sublimation radius 
but still outside the dust-free broad line region. A full discus- 
sion of this possibility, including realistic hot graphite grains 
spectra, is beyond the scope of this paper. 

Finally we note that the above suggested location, close or 
even inside the silicate dust sublimation radius, is also prob- 
lematic from the point of view of the derived covering factors. 
Our calculations assume no obscuration of one component by 
another yet such a location must obscure some of the radiation 
impinging on the torus. The location of these clouds imply 
that the radiation incident upon the torus should change and 
include part of the re-emitted radiation from the BB compo- 
nent. Given the unknown properties of the hot BB component, 
we cannot take this into account within the framework of the 
present work. 

5.4. Covering factors 

A major assumption of this work is that the entire MIR 
spectrum, after starburst subtraction, is reprocessed AGN ra- 
diation. This can be used to deduce the covering factor of 
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FIG. 10. — Distributions of different covering factors (arrows denote mean 
values). The total apparent CF is the the ratio between the integrated MIR 
luminosity and the bolometric luminosity. /((') is the ratio between the inte- 
grated, angle dependent, torus MIR luminosity and the bolometric luminosity 
of the source. This differs slightly from fh due to the anisotropic radiation of 
the torus. Note that PG 0157+001 and PG 0050+124 have been omitted since 
their CF exceed unity (see text). 



the central source by the three components (see also S08 and 
references therein; Maiolino et al. 2007). Regarding the BB 
and the NLR components, this factor is obtained directly from 
the comparison of their deduced luminosities with Lb i- The 
covering factor of the torus, fi, is different because of its non- 
isotropic radiation. For a given model, fi is calculated from 
the torus parameters. To demonstrate these differences, we 
compare in Table [6] fa, f(i), and the covering factors of the 
NLR and the BB components. The total covering factor can 
be defined in two ways: 1. Apparent total covering factor, 
defined as L b s (2 — 100) /Lbol- This is similar to the num- 
ber used in all earlier investigations of this type. Since the 
Spitzer-IRS spectra extend only to ~35 jJ.m, we use the mean 
"intrinsic" AGN SED of Netzer et al. (2007), obtained from 
the starburst subtracted QUEST spectra, to calculate the ratio 
between L bs(2 — 100) and L bs(2 — 35). This ratio is 1.072 
and thus we multiply L b s (2 — 35) of each object by this factor 
to get its apparent total covering factor. 2. Real total cover- 
ing factor, calculated by summing together fa and the NLR 
and BB covering factors. The various covering factor distri- 
butions are shown in Fig.fTOl 

The mean value of fa is 0.27. This value is lower than the 
/(/), for which the mean value is 0.34. The difference is due 
to the anisotropy of the torus radiation and the fact that, for the 
present sample of type-I sources, the torus inclination angle is 
small. The mean covering factor of the NLR component is 
0.07. This is in agreement with other estimates that are based 
on narrow emission line imaging and line equivalent width 
measurements, which suggest NLR covering factor < 10% 
(e.g., Baskin & Laor 2005; Netzer & Laor 1993). The mean 
covering factor of the BB component is 0.23. The mean real 
total covering factor in our sample is 0.57, which is slightly 
lower than the mean apparent total covering factor of 0.59. 

The dependence of the different covering factors on bolo- 
metric luminosity is seen in Fig. QTJ Looking at the entire 
sample, we find no significant correlation of bolometric lu- 
minosity with any of the covering factors. However, two of 
the sources, PG 1448+273 and PG 1244+026, seem to be af- 
fecting the entire correlation in a way that their removal from 
the sample makes the luminosity CF correlations significant. 
These sources are marked with different symbols on the dia- 
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FIG. 11. — Correlations between covering factors and bolometric luminosi- 
ties (left panels) and between torus mass and bolometric luminosity and the 
mass of the black hole (right panels). The square symbols represent sources 
with the lowest Mgn in the sample, see text for explanation 

gram. The reason why these objects should be omitted from 
the sample is not clear. We note, however, that these are 
the sources with the lowest Mbh thus may represent differ- 
ent properties. In particular, they are likely to be the highest 
amplitude variables and the highest accretion rate BHs thus 
their derived Lbol mav be wrong. Visual inspection of the bot- 
tom left panel shows that, without the two sources, the geo- 
metrical covering factor (fa+BB) decreases with Lbol and its 
highest value, at the low luminosity end, is roughly 0.7, in 
agreement with typical assumed distributions between type-I 
and type-II AGN in the local universe. 

Several earlier studies found clear anti-correlation between 
covering factor and Lb i (e.g., Wang et al. 2005; Richards et 
al. 2006; Maiolino et al. 2007; Treister et al. 2008). The 
decrease in covering factor translates, within the unification 
scheme, to a decreasing fraction of obscured AGNs as a func- 
tion of luminosity. All earlier studies were based on the total 
MIR emission and took no account of the various contributers 
to the flux and the fact that the geometrical covering factor 
cannot be obtained, directly, from L(IR)/Lboi- Our results are 
consistent with these studies but take into account, more ac- 
curately, the exact torus geometry and the needed translation 
between observed IR flux and geometrical structure. They 
show that ~ 7% of the derived covering factor is due to NLR 
clouds and properly account for the inclination angle distri- 
bution within the clumpy torus scenario. Taking the values 
found here we conclude that for Lb i— 5 x 10 45 ergs _1 there 
is roughly 1 : 1 ratio between type-I and type-II sources (note 
that the NLR contribution is subtracted from the total covering 
factor). Both Maiolino et al. (2007) and Treister et al. (2008) 
find a somewhat larger value for this bolometric luminosity. 
However, estimates of this ratio based on X-ray surveys are 
similar to the one found here and consistent across the entire 
luminosity range in our sample (e.g., Hasinger et al. 2004; 
Treister & Un-y 2006). 

The physical mechanism responsible for the decrease of 
covering factor with Lbol is still undetermined. One possi- 
bility is a receding torus mentioned above, where higher lu- 
minosity implies larger dust sublimation distance, and hence 
an obscuring structure that is located farther away from the 
center. In this scenario, the obscuring structure must have a 
constant height. In the clumpy torus model this scenario cor- 
responds to a decrease of a with Lb i- Indeed, we find that a 
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TABLE 6 

Covering Factors 



Object 


Real total CF 


Apparent total CF 


Real torus CF (/ 2 ) 


m 


NLRCF 


BB CF 


PG 2349-014 


0.70 


0.67 


0.30 


0.38 


0.09 


0.30 


FCj 2251+1 13 


0.30 


0.28 


0.09 


0.08 


0.05 


0.16 


FCj 2214+139 


0.36 


0.34 


0.13 


0.15 


0.03 


0.21 


PCj 1 700+518 


0.58 


0.57 


0.37 


0.42 


0.05 


0.17 


FCj 1626+554 


0.30 


0.33 


0.18 


0.21 


0.00 


0.11 


FCj 1617+175 


0.72 


0.73 


0.36 


0.44 


0.03 


0.34 


FCj 1613+658 


0.81 


0.78 


0.45 


0.60 


0.02 


0.34 


FCj 1448+273 


0.24 


0.24 


0.15 


0.16 


0.03 


0.06 


FCj 1440+356 


0.74 


0.79 


0.36 


0.45 


0.04 


0.35 


FCj 1435-067 


0.55 


0.54 


0.32 


0.34 


0.07 


0.17 


FCj 1426+015 


0.68 


0.69 


0.33 


0.46 


0.06 


0.29 


Di^ 1/111 1/1/11 

FCj 1411 +442 


0.92 


0.95 


0.48 


0.59 


0.03 


A A 1 
0.41 


PCj 1309+355 


0.40 


0.42 


0.16 


0.20 


0.10 


0.14 


FCj 1302-102 


0.77 


0.81 


0.40 


0.56 


0.14 


0.22 


pr. 1 9444-096 


44 


49 


U.JJ 


n 46 

U.'+O 


u.uo 


o 


PG 1229+204 


0.74 


0.86 


0.31 


0.41 


0.13 


0.29 


PG 1126-041 


0.91 


1.00 


0.37 


0.53 


0.16 


0.37 


PG 1116+215 


0.65 


0.63 


0.23 


0.28 


0.04 


0.38 


PG 1004+130 


0.33 


0.32 


0.20 


0.24 


0.05 


0.08 


PG 1001+054 


0.58 


0.64 


0.16 


0.21 


0.12 


0.30 


PG 0953+414 


0.51 


0.51 


0.10 


0.12 


0.06 


0.35 


PG 0838+770 


0.77 


0.80 


0.48 


0.65 


0.09 


0.20 


PG 0026+129 


0.45 


0.45 


0.14 


0.16 


0.07 


0.24 


PG 0157+001 


1.05 


1.08 


0.53 


0.80 


0.42 


0.11 


PG 0050+124 


0.89 


1.09 


0.63 


0.90 


0.03 


0.23 


B2 2201+31 A 


0.20 


0.22 


0.09 


0.12 


0.02 


0.09 


Mean 


0.57 


0.59 


0.27 


0.34 


0.07 


0.23 



is anti-correlated with the bolometric luminosity. 

We also find a positive correlation between the inclination 
angle of the torus and the bolometric luminosity. This is also 
consistent with a receding torus. An obscuring structure lo- 
cated at a larger distance, in this scenario, implies a larger 
solid angle through which an unobscured view to the center 
is possible. Thus, the luminosity of the source determines the 
possible range of inclination angles in which the model would 
be consistent with a type-I AGN SED. Since there is no pref- 
erence to a specific angle within that range, the correlation is 
due to the expansion of the possible range in i with increasing 
luminosity. 

5.5. Torus mass 

Given the torus parameters, we can estimate the torus mass 
from 



A/ """ 1 \0\m{a)N N H ^{^) r/,(F.. (9, 



Ma 



where /Vh.23 is the column density of a single cloud in units 
of 10 23 cm obtained from Ty and the assumed dust-to-gas 
ratio. I q is a function of Y and has the values of 1, Y /(2lnY) 
and |F for q = 1,2, and 0, respectively. 

The torus mass shows a clear correlation with both the bolo- 
metric luminosity and the mass of the central black hole (see 
Fig.fTTI). These correlations are expected since the bolometric 
luminosity appears in the torus mass calculation and, in gen- 
eral, Lbol ^A^bh- No correlation is found between the torus 
mass and L/L^. We also checked for correlation between 
Moms /Lbol and Lbol and found none. 

The mass of the torus ranges from 8.45 x 10 4 to 3.12 x 
1O 7 M for the silicate dust sublimation radius. The ratio be- 
tween the torus and the BH masses is ranges from 4.5 x 10~ 4 



to 0.08 and does not show any correlation with the source lu- 
minosity. 

Assuming a mass-independent BH accretion efficiency of 
tj = 0.1, we can calculate the time it would take for the torus 
mass to be completely accreted by the BH. This is an ex- 
tremely short duration of about ~ 10 s yr. This time scale 
is similar for all sources and is independent of torus mass, 
source luminosity and any of the other properties. The typical 
e-folding time for a BH mass of 10 8 M Q and the above T] is 
about 4 x 10 7 (L/LEdd) 'yr- This mean that, if all accretion 
to the center goes through the torus, there are more than 100 
short accretion episodes per one e-folding growth time. 

The above result indicates at least two different possibil- 
ities: 1. A constant replenishment of the torus material by 
larger scale mass inflow through the plane of the torus. Such 
a scenario, requires a constantly changing torus structure. 2. 
Matter accreting onto the BH does not originate in the torus. 
Such a scenario has been suggested by Elitzur & Shlosman 
(2006), who assumed that mass from the host galaxy arrives 
directly into the accretion disk. Consequently, wind from the 
accretion disk forms the toroidal obscuration. As long as the 
large scale mass infall is sufficiently large and continuous, the 
wind sustains a steady outflowing clumpy structure, which is 
the clumpy torus discussed in this work. In this case, the small 
mass of the torus is not directly related to the growth time of 
the BH. 

We are grateful to Dieter Lutz, Mario Schweitzer, Benny 
Trakhtenbrot and Ido Finkelman for useful discussions. We 
also thank Benny Trakhtenbrot for retrieving the SDSS data 
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this work has been provided by the Israel Science Foundation 
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